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ABSTRACT 
In this paper we discuss a dynamic flux-transport dynamo model that includes the feedback of 
the induced magnetic field on differential rotation and meridional flow. We consider two different 
approaches for the feedback: mean field Lorentz force and quenching of transport coefficients such as 
turbulent viscosity and heat conductivity. We find that even strong feedback on the meridional flow 
does not change the character of the flux-transport dynamo significantly; however it leads to a signif- 
icant reduction of differential rotation. To a large degree independent from the dynamo parameters, 
the saturation takes place when the toroidal field at the base of the convection zone reaches between 
1.2 an 1.5 T, the energy converted into magnetic energy corresponds to about 0.1% to 0.2% of the solar 
luminosity. The torsional oscillations produced through Lorentz force feedback on differential rotation 
show a dominant poleward propagating branch with the correct phase relation to the magnetic cycle. 
We show that incorporating enhanced surface cooling of the active region belt (as proposed by Spruit) 
leads to an equatorward propagating branch in good agreement with observations. 
Subject headings: Sun: interior — rotation — magnetic field — dynamo 



1. INTRODUCTION 



Flux-transport dynamos have proven to be success- 
ful for modeling the evolution of the large scale so- 
lar magne tic field IWang fc Sheelev"1991^, 'Durnev"199a 
^houdhuri '^t all [l995: Dikoati & Charbonncau 19q1 
Kiiker ct al. 2001: Dikoati & Gilman 2GGlt lDikpati et all 
2004: Dikpati 2005: Dikoati et al. 200(1 . In a fiux- 



transport dynamo the equatorward propagation of the 
magnetic activity belt (butterfiy diagram) is a conse- 
quence of the equatorward transport of magnetic field at 
the base of the convection zone by the meridional fiow. 

However, all studies so far have addressed the trans- 
port of magnetic field by the meridional circulation in a 
purely kinematic regime. The toroidal field strength at 
the base of the solar convecti on zone inferred from studies 
of rising magneti c flux tubes ( Choudhuri & GilmanT987; 
[ Fan et alJ 119931 iSchiissler et al. 1994: Cahgari et a l. 
II 995111 9981) is around 10 T (100 kC) and thus orders of 
magnitude larger than the equipartition field strength es- 
timated from a meridional fiow velocity of a few ms~^. 
Therefore it is crucial for fiux-transport dynamos to in- 
clude the feedback of the Lorentz force on the meridional 
fiow. 

In order to be able to address this question it is nec- 
essary to incorporate a model for the solar differential 
rotation and meridional flow into a dynamo model and 
allow for the feedback of the Lorentz force on differential 
rotation and meridional flow. Differential rotation and 
meridional flow have been addressed in the past mainly 
through two approaches: 3D full spherical shell sim- 
ulations fGlatzmaier & Gilman" "1982"; "Gilman & Miller! 
Ease; Miesch et al. 2000; Brun fc Toomre 2002) and ax- 
isymmetric mean field models ( Kitchatinov fc Riidiger 
Il993l IT995t IRiidiger et aP IT99a .Kiiker fc Stix. .200 U . 
While the 3D simulations have trouble reproducing a 
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consistent large scale meridional fiow pattern (pole- 
ward in the upper half of the c onvection zone) , as 
it is inferr ed by hclioseismology (Brau n fc FanI Il998t 
iHaber et al. 2002: Zhao fc Kosovichev 2004J) . such a flow 
is a common feature in most of the mean field models. 

In this paper we build upon the differential rotation 
model presented in RemDcl (2005b) and combine it with 
the axisymmetric mean field dynamo equations. The 
magnetic field is allowed to feed back on differential rota- 
tion and meridional flow through the mean field Lorentz 
force and quenching of turbulent viscosity and heat con- 
duct ivity. We have addressed already in a previous pa- 
per (Rempel 2006) the feedback of magnetic field on 
the meridional flow by imposing a static toroidal mag- 
netic and including magnetic tension and quenching of 
turbulent viscosity in the differential rotation model. 
We found that a signiflcant feedback can be expected 
if the toroidal magnetic field strength is around 3 T or 
above equipartition. In this paper we do not impose a 
toroidal field, but rather solve the induction equation 
to ob tain a time de pendent magnetic field. Going be- 
yond iRempell l)2006j) , we also include the feedback on 



differential rotation leading to a solar cycle variation 
of the ro tation rate, which is k nown a s torsional os- 
cillations llHoward fc LabontdlT98 0: Too mre et al.ll2000l: 
[Howe et abl I200nt lAntia fc Ra,suri2001; .Vorontsov et alJ 

i2oo2HHowe ^tmmmi . 

A very simi lar approach was taken before by 
iBrandenburg et al. (1990, 1991, 1992). In their model 
they were solving a mean field differential rotation model 
parallel to the dynamo equation to obtain a 'dynamic' 
dynamo. The main difference in our approach is that we 
focus on flux-transport dynamos, whereas their work de- 
scribed mainly afi-dynamos with the advection of field 
by a meridional flow playing only a secondary role. Given 
recent developments in solar dynamo theory, showing 
that flux transport dynamos are very s uccessful in rep ro- 
ducing most of the observed features jDikpatil l2005f) it 
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is important to evaluate to which degree these dynamos 
change operation if dynamic feedback is considered. 

Lorentz force feedback has been considered in mean 
fi eld dynamo models in various levels of sophistication: 

'Kii ker e t al." (1999) used a model including both macro 
(mean field Lorentz force) and micro (quenching of A- 
effect) feedback to evaluate to which extent grand min- 
ima can be produced through feedback on differential 
r otation. 

'Covas et al." ("200^, "^OOl", "2005^ considered feedback on 
differential rotation in a classical afJ-dynamo model solv- 
ing a simplified momentum equation including the mean 
field Lorentz force and a diffusive relaxation term for 
the longitudinal flow velocity perturbation. They were 
able to reproduce the basic features of the observed so- 
lar torsional oscillation pattern (equatorward and pole- 
ward propagatin g branch). A similar approach has bee n 
taken befor e bv [Moss fc Brook^ l)200o|) : iTobiaJ l)1996|) : 
lYoshimural (|1981f l. ^^^_^^^ ^_^ 

In contrast to this the models of Schiissler ("197^; 
[Brandenburg et al . (1990, 1991, 1992)- Jennings (,199J; 
IMoss et all l|1995|) : Muhh et al. ( 199j) incorporate the 
full momentum equation, allowing also for magnetically 
driven meridional motions. The main focus of their work 
was on understanding the non-linear saturation of the 
dynamo. 

The consideration of the macroscopic mean field 
Lorentz force, common for all models listed above, is 
also known in the litera ture as 'Malkus-Proctor-effect' 
l|Malkus fc ProctoillT97l . 

A different ap pr oach has been used by 
iKitchatinov fc PipinI l)1998j) and iKitchatinov et all 
l)1999D . who considered feedback through quenching of 
the A-effect (turbulent angular momentum transport 
driving differential rotation). Their work focused on 
understanding torsional oscillations as well as the 
possibility of producing grand activity cycles through 
this type of feedback. 

This paper is organized as follows: In section [3 we ex- 
plain the physics included in the non-kinematic dynamo 
model. Section 13 shows the results of the non- kinematic 
dynamo runs including a detailed analysis of the energy 
flows within the model. Section 0] focuses on the proper- 
ties of the torsional oscillations produced by the model 
and compares them to results obtained by helioseismol- 
ogy. Section [S] discusses the choices we make for various 
parameters of the mean field model and their impact on 
the solutions presented here. In section IHl we summarize 
the main results of this investigation and discuss them 
in the context of solar dynamo models. 

2. MODEL 

In this paper we utilize the mean field dif ferential 
rotation and meridional circulation model of iRempell 
l)2005b() and couple it with the axisymmetric mean field 
dynamo equations. The computed differential rotation 
and meridional flow are used to advance the magnetic 
field by using a Babcock-Leighton flux transport dynamo 
model. The computed magnetic field is allowed to feed 
back on differential rotation and meridional flow through 
the Lorentz force and quenching of turbulent viscosity 
and thermal heat conductivity. We are solving the ax- 
isymmetric MHD equations including parameterizations 
of processes on the (unresolved) convective scale (mean 



field approach). We introduce in the induction equation 
the vector potential for the poloidal field to satisfy the 
constraint 'V ■ B — 0. 



dgi 1 d 



1 d 



[suveveQo) (1) 



ot r'^ or ^ ' r sm oO 

dVr dVr V0 dVr , Wg 3 Ptot , Pmag , Si 

'TTT^"^^^ 1^+ 1^ + -9+—9 

ot or r OO r or qq 7po 7 

+ (2r!or!i + nl) r sin2 6+ — (F^ + F^) (2) 

Qo 

dvg__^ dvg Vg dvg VrVQ Id Ptot 

dt ^ dr r 89 r r 89 qq 

+ (2f7or2i + nj)rsm9cos9+ — (F^ + F^) (3) 



8VI1 Vr 8 

8t r^ 8r 

Vg 8 



■siii^ 9 89 



[r^ino + Qi)] 

[sin^ 9{flo + ni 



1 



Qor sin 9 



{Fl + Fi 



(4) 



dsi dsi Vg 8si 7*5 , 7- 1^ 

8t 8r r 89 Hp Pq 

■ ^ V-{KtQoTaVsi)^^^ijt{V^Bf (5) 



QoTq 



Po 



8B^ 18, ^ 1 9 , , ^ dfli 

— — = — [rVrB^) — {vgB^) + rsm9Br-^— 

8t r 8r r 89 8r 

+ sin.i?.— +,.^A-^---^ji?, 

4-i§^^(.i?.) + l§-i^|(sm.i?,) (6) 
r 8r 8r r^ 89 sm 9 89 

dA Vr 8 vg d 

8t r 8r r sm 9 89 



+Vt[A 



1 



(r sin ( 



A + S{r,9,B^) 



(7) 



Here go and po denote the (spherical symmetric) ref- 
erence state stratification, Hp = po/{gog) the pressure 
scale height, gi and pi perturbations around the ref- 
erence state caused by differential rotation and merid- 
ional flow. Since these perturbations are small com- 
pared to the reference state values, the equations are 
linearized assuming gi <C go and pi <C po- ^0 denotes 
the rotation rate of the core, Oi the differential rota- 
tion with respect to the core in the convection zone. 
The quantity si = pi/po — 7Pi/po denotes the dimen- 
sionless entropy perturbation (normalized by the heat 
capacity Cy). For the reference state we use an adia- 
batic polytrope assuming a gravity varying ~ r~^; how- 
ever, small perturbations from adiabaticity are consid- 
ered in the entropy equation through the third term 
~ (5 = V — Vad- The quantity pniag denotes the magnetic 
pressure, ptot — Pi + Pmag the total pressure. The buoy- 
ancy term in Eq. (0) has been written in a way to sep- 
arate the magnetic buoyancy from non-magnetic buoy- 
ancy, assuming |V — Vad| <^ 1- Since the most unstable 
modes driven by magnetic buoyancy are typically non- 
axisymmetric, the description of magnetic buoyancy in 
our axisymmetric model is not necessarily very realistic. 
Therefore we will discuss later also simulations with mag- 
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netic buoyancy switched off by ignoring the term ~ Pmag 
in Eq. @. 

The quantity F = l//ioV • (BB) denotes magnetic 
tension, where the poloidal magnetic field foUows from 
the vector potential A used in Eq. (O through 
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For computing the Lorentz force we consider here only 
the magnetic mean field contribution. Formally an addi- 
tional contribution to the stress tensor ~ {B'B') caused 
by the turbulent magnetic field exists; however the con- 
tribution of these terms in detail is not well under- 
stood. We emphasize that considering only the mean 
field Lorentz force leads to a model that is energetically 
consistent in the way that the energy extracted from 
differential rotation and meridional flow is flowing into 
the reservoir of magnetic energy via the induction equa- 
tion. Considering additional Lorentz force terms in the 
momentum equation requires that these terms only re- 
distribute momentum (otherwise additional terms in the 
induction equation would be required as well to be en- 
ergetically consistent). Thi s effect referred to as visco- 
elasticity by some authors (iLonecope et alJl2003[l could 
parametrized to some extent as an additional source of 
viscous stress. Another possible feedback of these terms 
is quenching of turbulent transport processes such as tur- 
bulent angular momentum transport, viscosity and heat 
conductivity, which we will consider later in the form: 
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The quantity F'^ denotes turbulent viscous stresses, 
including turbulent transport of angular momentum (A- 
effect), Q considers the associated viscous dissipation in 
the entropy equation. In our model the A-effect is the pri- 
mary driver of differential rotation and meridional flow, 
while the profile of the differential rotation (the devia- 
tion from the Taylor-Proudman state) is a consequence 
of a latitudinal entropy gradient. The latitudinal entropy 
gradient follows in our model self consistently from the 
inclusion of a subadiabatic tachocline. For further de- 
ta ils conce r ning th e differential rotation model we refer 
to'RemiDcl' ('2005b') and section [Q 

Eqs. ^ and Q are the axisymmetric mean field dy- 
namo equation, including transport of magnetic field by 
meridional flow, shear by differential rotation, magnetic 
diffusion and induction by a Babcock-Leighton surface 
a-effect, parametrized through the poloidal source term 
S (r, 0, i?$). We will discuss the dynamo model in more 
detail in section IT^ 

For the solar parameter range the system defined by 
Eqs. Q) - (0 is in the regime of highly subsonic flows, 
which introduces a significant CFL time step constraint 
if solved explicitly. We do not use here the anelastic 
approximation as many others do; instead we use an ap- 
proach of artificially reducing the speed of sound to val- 
ues that do not impose a significant time step constraint 



(compared to the diffusive time step) but still ensure that 
the flows remain highly subsonic. Given the fact that the 
Mach number of the meridional flow in the bulk of the 
convection zone is around 10~^ a decrease of the speed 
of sound by a factor ~ 100 is found to have no impact 
on the solution and is used for most of the results shown. 
Formally this is achieved by changing the equation of 
continuity to 
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with C ~ 0.01, which reduces the speed of sound by a fac- 
tor (. This approach is equivalent to increasing the base 
rotation rate and scaling up all other variables to main- 
tain the prope r relations betw een the different terms in 
the equations l)R,emDell|2005bD . The modified equations 
are solved using a MacCormack scheme. We emphasize 
that this approach is only feasible for the axisymmetric 
system, where the much faster speed of rotation does not 
enter the CFL condition. 

2.1. Differential rotation, meridional flow reference 

model 

In this paper we use a reference model t hat is very close 
to model 1 discussed in iRempell l|2005bj) . We made the 
following changes in order to obtain a meridional fiow 
pattern that leads to a solar like dynamo period and 
a confinement of magnetic activity close to the equator: 
We use a value of the parameter n defining the latitudinal 
profile of the A-effect of 3, an amplitude of the A-effect 
of Aq — 1, and a value of turbulent viscosity and heat 
conductivity of 3 x 10®m^s~^. Aq determines primarily 
the amplitude of differential rotation, while a change of 
vt and Kt (keeping vt/ K,t constant) adjusts the meridional 
flow speed. 

Fig. n] a,b) shows the differential rotation and Fig. ^ 
c,d) the meridional flow of the reference model. Note that 
in our model most of the radial shear is located beneath 
the base of the convection zone rbc — 0.71 i?©. This is 
due to the fact that the differential rotation is driven by 
the A-effcct within the convection and uniform rotation 
is imposed at the lower boundary at r = 0.65 i?©. As 
a consequence a viscous shear layer forms between both 
regions, which has the largest shear rate where the tur- 
bulent viscosity is assumed to be small (below rbc). Also 
the meridional flow does not show a significant penetra- 
tion below The , due to the subadiabatic stratifica tion and 
the sig nificant drop of turbulent viscosity (see IRempell 



the sig r 
Il2005bl) 



for a detailed discussion). 



2.2. 



Flux-transport dynamo model 

We use a flux-transport dynamo rnodel si milar to the 
approach of Dikpati & Charbonneau (1999), except for 
the fact that differential rotation and meridional flow 
are not prescribed, but computed by the model de- 
scribed above. The additional dynamo parameters we 
have to specify are the profile of the magnetic diffusivity 
rjt and the functional form of the poloidal source term 
S (r, 9, _B$) in Eq. [3 The turbulent magnetic diffusivity 
rjt is specified as function in radius, given by 

■r]t^'ric + fc{r) [rjhc ~ Vc + fczir) (r/cz - ??bc)] (13) 
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Fig. 1. — Differential rotation and meridional flow of the reference model, a) Contours of f2, b) radial profiles of f2 for the latitudes (top 
to bottom) 0°, 15°, 30°, 45°, 60°, and 90°. c) shows the stream function of the meridional flow, d) the flow profile at 30° latitude. The 
dotted line indicates the base of the convection zone at r^^ = 0.71 Rq. Note that in our model most of the radial shear is located beneath 
rj-,;,, while the meridional flow shows only very little penetration beneath r^,.. 




Fig. 2. — Profile of rjt used for the dynamo simulation, rjt drops 
by two orders of magnitude within the convection starting from a 
surface value of lO^m^s"^. There is an additional drop by one 
order of magnitude beneath the base of the convection zone. 




Fig. 3. — Butterfly diagram of toroidal magnetic field (contour 
lines of S#,bc) and radial fleld close to surface r = 0.985 Rq (color 
shades). The maximum field strength of _B$ is 1.28 T, the maxi- 
mum field strength of Br is 0.01 T. 
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The function fcz{r) determines the profile within the 
convection zone, while fdr) ensures a significant drop 
of rjt beneath rtc- ?7c determines the core diffusivity, ?7bc 
the diffusivity at the base of the convection zone, and r/cz 
the diffusivity in the upper half of the convection zone. 
For this discussion we use a profile with the parameters 



77c = lU"m*-s % rjhc = lU"m's , rjcz 

rcz = 0.875 i?0, dcz = 0.05 i?0, rbc = 0.71 i?0, dbc = 

0.0125 i?0. The profile is shown in Fig. [3 We use here 

for reasons of numerical stability a value of rjc that is 

significantly larger than molecular resistivity. However, 

the influence of rye on the solution is found to be very 

weak. 

Our model uses a non-local B abcock-Leighton a-effect 
('Babcock""1 961t lLeisf hton""1969") in which the source term 
S {7-,6, B<i,) at the surface is dependent on the toroidal 
field strength at the base of the convection_ zone aver- 
aged over the interval [0.71 i?©, 0.76 i?©], i?$,bc- The 
functional form of S is then 



S (r, 0, B$) = ao B^M^) /^(r) 5„(0) , (16) 
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(17) 
(18) 
(19) 

We use here d^ = 0.05 i?©, which confines the poloidal 
source term above r = 0.935 i?0, peaking at rmax- The 
function h{r) is an averaging kernel with J """^ Ar h{r) — 
1. We use a parabolic profile vanishing at 0.71 Rq and 
0.76 i?0 with peak at 0.735 i?©. 

The boundary condition is A = 5$ = at the pole and 
dA/d0 = Bg, = at the equator, which selects the dipole 
symmetry for the solution. _B$ vanishes at both radial 
boundaries. A vanishes at the inner boundary and the 
poloidal field is assumed to be radial at the top boundary. 
The dynamo part of o ur code has been compared inte n- 
sively with the code of lDikpati fc Charbonne;iul l)1999(l . 

We emphasize that the main goal of this paper is a 
fundamental understanding of dynamical effects caused 
by the feedback of the dynamo generated field on dif- 
ferential rotation and meridional flow, rather than a de- 
tailed model of the solar dynamo. We have chosen the 
dynamo parameters and the parameters of differential ro- 
tation model such that the solutions show a reasonable, 
but not detailed agreement in terms of butterfly diagram 
and dynamo period as well as amplitude of differential 
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rotation and meridional flow with observations. We re- 
strict our simulations to one hemisphere and impose the 
dipole symmetry through our equatorial boundary con- 
dition. We do not try to address the parity issue in this 
m odel, which has bee n done for flux-transport dynamos 
bv lDiknati fc Gilmanl pOOl) . 

Before we discuss dynamic solutions, we present here 
as reference a kinematic dynamo solution computed with 
the differential rotation and meridional flow shown in 
Fig. n We use for the a-effect and amplitude of ao = 
0.125 ms^^ and include a-quenching with a quenching 
field strength of IT (10 kG). The butterfly diagram 
computed from the averaged toroidal field (i?$,bc) and 
the radial field at the 'surface' r — 0.985 i?0 is shown 
in Fig. 121 The dynamo period for this setup is around 
19 years, the maximum toroidal field at r = 0.735 i?0 is 
1.28 T, the maximum field strength of the radial field at 
r = 0.985 i?0 is 0.01 T. 

The butterfly diagram shows equatorward propagating 
activity belts starting around 50° latitude and having 
their peak field strength at around 40° latitude. The po- 
lar reversal of the poloidal field takes place during max- 
imum activity with the toroidal field in low latitudes, 
changing sign from negative to positive while the low lat- 
itude toroidal field is positive as found in observations. 

3. DYNAMO MODEL WITH LORENTZ FORCE FEEDBACK 

3.1. General solution properties 

In this section we discuss results obtained with full 
Lorentz force feedback on differential rotation and merid- 
ional flow. Since Lorentz force feedback introduces 
enough non-linearity to saturate the dynamo, it is not 
necessary to include a-quenching as typically done in 
kinematic models. We present here three models vary- 
ing the value of the a-effect in order to show different 
regimes in terms of intensity of the Lorentz force feed- 
back. The ao values used are 0.125, 0.25, and 0.5 ms~^. 
We also present results computed from a model with a 
magnetic diffusivity reduced in the bulk of the convection 
zone from yycz = 10^ m^ s""'^ to jycz = 5 x 10^ m^ s~^ . 

Figure 21 shows the model with a^ = 0.125 ms^^. The 
top panel displays the radial magnetic field close to the 
surface (r = 0.985 i?0), the middle panel the meridional 
flow at r = 0.735 i?0, and the bottom panel the tor- 
sional oscillations at r = 0.985 i?©. In all three panels 
the toroidal fleld contours at r = 0.735 i?0 are indicated 
(solid: positive values, dashed: negative values). Com- 
pared to Figure|31the butterfly diagram shows major dis- 
tortions caused by magnetic buoyancy breaking up the 
layer of toroidal field at the base of the convection zone. 
We present a more detailed discussion of the role of mag- 
netic buoyancy in subsection l3.4l The radial surface field 
is close to that in Figure 13 The equatorward meridional 
flow at r = 0.735 Rq shows a variation of around 30% of 
the mean flow amplitude (around 2ms~"'^) in anti-phase 
with the toroidal field intensity, caused by the influence 
of the magnetic tension of the toroidal field. This feed- 
back is not strong enough to switch off the equatorward 
transport of magnetic field and therefore does not influ- 
ence the flux transport dynamo significantly. This is in 
agreement with RemDcl (2006), who found that equator- 
ward transport of toroidal field is possible up to around 
3 T. The bottom panel of Figure 0] shows the tempo- 
ral variation of the differential rotation caused by the 




Fig. 4. — Dynamo solution with Lorentz force feedback and no 
o-quenching. All dynamo parameters are the same as for the ref- 
erence solution (with a-quenching) shown in Fig. J^J- Top panel: 
Butterfly diagram of toroidal magnetic field (contour lines) and 
radial field close to surface r = 0.985 -R0 (color shades). The max- 
imum toroidal field strength is 1.17 T, the maximum radial field 
strength at the surface is 0.014 T. Middle panel: Meridional fiow 
{vg, color shades) at r = 0.75 -R0. The maximum flow velocity is 
2.27 ms~^. Bottom Panel: Torsional oscillation pattern (fi — H) at 
r = 0.985 Rq. The maximum amplitude is 1% of the core rotation 
rate corresponding to a variation of about 4 nHz. The torsional os- 
cillation pattern at r = 0.735 Rq is very similar, but has a slightly 
reduced amplitude. 



Lorentz force feedback, also known as torsional oscilla- 
tions. We find in our model mainly a poleward propa- 
gating oscillation pattern, starting at mid latitudes. The 
amplitude close to the surface is around 1% of the core 
rotation rate, which corresponds to roughly 4 nHz. We 
discuss torsional oscillations in detail in section 0] 

Fig. El shows 3 snapshots of the magnetic field (field 
lines of poloidal field and toroidal field strength as color 
shades), corresponding to i = 0, i = 3.25 and t — 6.5 
years in Fig. 01 The middle and bottom panel show 
the breakup of the magnetic layer caused by magnetic 
buoyancy force. 

Fig. shows the dynamo solution (butterfly diagram 
and radial field close to surface) for the cases with in- 
creasing value of uq. In a kinematic model with a- 
quenching as non-linearity, an increase of a increases the 
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Fig. 6. — Butterfly diagram (contour lines) and radial field at 
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in Lorentz force feedback leads to a confinement of magnetic ac- 
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magnetic buoyancy. 



field strength, but does not influence the shape of the 
solution in great detail (with increased quenching the 
average profile of a changes, which has a slight influ- 
ence of the solution). In our model the toroidal field 
shows only minor changes (1.2, 1.4, and 1.2 T for the 
cases with ao values of 0.125, 0.25, and 0.5 ms~^, re- 
spectively), however the profile of the solution changes 
significantly through the changes in meridional flow and 
differential rotation. On average stronger feedback leads 
to a concentration of magnetic activity to lower latitudes, 
which is the consequence of a strong cycle variability of 
the meridional flow. While the average flow is roughly 
the same in all cases, the meridional flow tends to be 
more concentrated toward the equator during the phase 
of the cycle when the poloidal fleld is transported down- 
ward, leading to the production of toroidal fleld at lower 
latitudes, too. Even though interesting, this feature has 
most likely no relevance to solar dynamos, since the de- 
gree of feedback required is not observed (e.g. torsional 
oscillations with an amplitude of around 20 nHz). When 
considering solutions with such a large degree of feedback 
mag neti c buoyancy also plays an important role (subsec- 
tion EHJ. 
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3.2. Energy flows within model 

In this section we discuss the dynamo solutions by an- 
alyzing the energy flows between the different energy 
reservoirs of the model. This allows us to understand 
the saturation mechanism of the non-linear dynamo on 
a more quantitative level. 

Figure [7| shows a schematic of the energy flows in our 
coupled dynamo-differential rotation model. For reasons 
of simplicity we consider here only the energy of the 
toroidal field, since the energy of the poloidal field plays 
only a minor role. The following equations describe the 
change of rotation energy _En, energy of meridional flow 
Emi and energy of toroidal magnetic fleld Eb'- 



dEn 
dt 

BEm 
dt 

OEb 
dt 



— Qk-Qy- Qc - 
= Qc — Qy - Qb 



iL 






■ Ql + Ql Q-q : 



where the energy reservoirs are given by 



1 



En= I dV^gos^n'' 



E 



M '■ 



dVgo- 



EB=fdV^ 
J 2^0 

and the exchange terms are given by 



g/ 



Q^ 



Qc 






dVvtgos 



dV s g^v,-, 



on 

dr 



fidn\ 
[r'de) 



.n^ 



(20) 
(21) 
(22) 

(23) 
(24) 
(25) 

(26) 

(27) 
(28) 



Q2 = - I dVnBp-V{sB^) 

Ql, — / dV — [RrrErr + 2RreEre 



+R<p<f,Eij,ij,] 



(29) 



(30) 



Qb = - dVv, 


^Qog— 

7 


(31) 


QM^JdV^Vrr^-VisB^) (32) 


Q.-J^v'^ 


\ dr ) ^ \r do ) 


(33) 



Here s = r sin 9 denotes the distance to the axis of rota- 
tion and 



dy = 47r 



dr 



-it/2 



dOr'^ sin 9 



(34) 



denotes the integral over the entire volume of the sphere 
from r = Tmin to r = Tmax- Wc emphasize that we 
solve our model only for the northern hemisphere but 
we compute from that the energy conversion for the en- 
tire sphere. The quantity v^ = {vr,vg,Q) denotes the 
meridional flow and Bp — {Br,Bg,0) the poloidal mag- 

n etic field . 

iRempell l^2005a^ showed already the derivation for the 
non-magnetic part of the system Eqs. H20|l - (|22(l . The 
exchange terms Eqs. (j^ - (j^ are not unique expres- 
sions since they can be transformed in the form: 



expressioni = V • flux + expression2 



(35) 



provided that the volume integral over the flux diver- 
gence vanishes. We use the appropriate closed bound- 
ary conditions for all variables except for _B$, however, 
the resistive flux across the upper boundary turns out 
to be negligible. We also emphasize that we defined a 
few of the excha nge terms with an opposite sign than in 
IRempell l|2005aD . Here all Q-- are positive so that the 
signs in Eqs. H2()(l - (|22(l clearly state which terms are 
sources and which are sinks for the corresponding energy 
reservoir. 

The third column of tablensummarizes the energy flow 
of the differential rotation reference model we use for all 
dynamo simulations. In this model the A-effect converts 
1.4% of the solar energy flux into rotation energy. 57.4% 
of this energy is lost directly through viscous dissipation 
of the differential rotation, while 42.5% is flowing into 
the reservoir of meridional flow by means of the Coriolis 
force. The major fraction of this amount returns into the 
reservoir of internal energy through work against buoy- 
ancy force, only a small fraction < 1% through viscous 
dissipation of the meridional flow. In the reference model 
the pole-equator difference in rotation rate is 27% of the 
core rotation rate. 

Columns 4 to 10 of table ^ summarize the results ob- 
tained from the dynamo models with different a and 
77cz values. The top portion of Table H shows the pole- 
equator difference in rotation rate, the amplitude of tor- 
sional oscillations, the maximum toroidal and radial field 
strength, average magnetic energy and fluctuation of 
magnetic energy over a dynamo cycle as well as max- 
imum magnetic energy and flux of one polarity at the 
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TABLE 1 
Summary of results 



quantity 



unit 



reference ao = 0.125nis 



ao = 0.25 ms 



ao 



0.5ms 



rjcz 



(f^eq - fipolc) /^O 




0.27 


0.21 


(0.21) 


0.15 


(0.15) 


0.1 


(0.12) 


0.15 


max(C — n)gQO 


[nHz] 






4.7 


(5.8) 


11.5 


(11.1) 


17.7 


(24.2) 


12.9 


max(Q - n)gQO 


[nHz] 






3.5 


(2.9) 


6.9 


(5.8) 


8.3 


(12.1) 


9.7 


max(i?$) 


[T] 






1.2 


(1.2) 


1.4 


(1.1) 


1.2 


(1.2) 


1.5 


max (Br) 


[T] 






0.014 


(0.015) 


0.022 


(0.021) 


0.027 


(0.047) 


0.02 


Eb 


[1031 J] 






2.8 


(3.5) 


4.6 


(4.4) 


5.1 


(4.6) 


4.9 


max[(£;s -Eb) /Eb] 








0.12 


(0.15) 


0.22 


(0.26) 


0.26 


(0.55) 


0.21 


max{£;+)b^ 


[1031 J] 






1.7 


(2.1) 


2.2 


(2.2) 


2.1 


(2.7) 


2.4 


max($+)bc 


[IQieWb] 






1.1 


(1.3) 


1.2 


(1.3) 


1.2 


(1.6) 


1.2 


Qa 


[Ee] 


0.014 


0.013 


(0.013) 


0.011 


(0.01) 


0.008 


(0.008) 


0.012 


Q? 


[Qa] 


0.574 


0.5 


(0.484) 


0.405 


(0.415) 


0.304 


(0.357) 


0.459 


Qc 


[Qa] 


0.425 


0.43 


(0.429) 


0.445 


(0.423) 


0.465 


(0.412) 


0.441 


Qf 


[Qa] 




0.069 


(0.087) 


0.149 


(0.162) 


0.232 


(0.231) 


0.1 


Qb 


[Qa] 


0.419 


0.41 


(0.404) 


0.414 


(0.384) 


0.423 


(0.357) 


0.409 


Qt' 


[Qa] 


0.005 


0.005 


(0.005) 


0.007 


(0.007) 


0.01 


(0.01) 


0.008 


Qi' 


[Qa] 




0.014 


(0.019) 


0.024 


(0.032) 


0.033 


(0.043) 


0.023 


Qv 


[Qa] 




0.083 


(0.104) 


0.172 


(0.192) 


0.268 


(0.272) 


0.123 


Qr, 


[^q] 






0.0011 


(0.0014) 


0.0019 


(0.0019) 


0.0021 


(0.0022) 


0.0015 



Note. — Summary of dyna mo s im ulat ions discussed. The top portion shows a few global properties of the solution, the bottom portion the enei 
transfer terms defined by Eqs. I26i - I33i . Values in brackets refer to models in which the magnetic buoyancy term is switched off (see subsection | 
for further details). The maximum of H — H and Br is evaluated at 0.985 iJ©, the maximum of B^ at 0.735 iJ©. max(i?J)bc and max($+)i,c are i 
maximum values of magnetic energy and magnetic flux integrated over the interval [0.71,0.76] Rq, considering only one toroidal field polarity. The va 
of Qa and Q,, (last row) are given relative to the solar luminosity, all the other energy exchange terms are relative to Qa. Note that the accuracy of i 
energy exchange terms is around 0.001, so the equilibrium relations Qa = Ql? + Qc + Ql' Qc = Q*^ + Qb + Qift and Q,, = Q^ + Q^^ are only fulfil 
within that error margin. 



base of the convection zone. The bottom portion shows 
the energy exchange terms Eqs. (|26|) to (|33|l averaged 
over a cycle (12 cycles in case of the irregular solutions 
shown in Fig. IH)). The quantity Qa is given in units of 
the solar energy flux and the following terms relative to 
Qa- The last line gives Q^ in units of the solar energy 
flux. 

As already mentioned in the previous section, increas- 
ing the value of a does not lead to a significant increase 
of the toroidal field strength. The same applies to the 
magnetic energy and the magnetic flux at the base of 
the convection zone. At the same time the equator-pole 
difference of the differential rotation decreases monoton- 
ically from 0.2filo in the case with ao = 0.125 ms^^ to 
0.1 r^o in the case with ao = 0.5 ms~^. For compari- 
son, the reference model has an equator-pole difference 
of 0.27ilo. This shows clearly that the saturation mecha- 
nism of the dynamo is the reduction of the shear through 
Lorentz force feedback on differential rotation. This be- 
comes also evident by looking at the energy exchange 
terms Q^ and Q^: While Q^ is decreasing and Q^ in- 
creasing, the sum of both relative to Qa remains roughly 
the same, meaning that in the dynamo solutions viscous 
stress is replaced by Maxwell stress. The energy transfer 
to meridional flow Qc shows only a weak increase with 
ao- The energy transferred in total into magnetic energy 
(and dissipated back to internal energy), Q^ increases 
relative to Qa from 8.3% to 26.1% with increasing ao. 
Even though we use a fixed parameterization for the A- 



effect in our model, the energy which is converted by the 
A-effect (Qa) decreases with increasing ao, since in Eq. 
(j^ also the shear enters. As a consequence, Q^ does 
not change that much in absolute units (last row in Ta- 
ble^) when ao is increased from 0.25 ms~^ to 0.5 ms~^. 
This is consistent with the almost constant value of the 
magnetic energy in both cases. 

A different way to look at the saturation is as follows. 
In the model with ao = 0.5 ms~^ the radial surface field 
(a measure for the poloidal field strength) almost dou- 
bled compared to the model with ao — 0.125 ms^^. At 
the same time equator-pole difference of fl is reduced to 
roughly half the value, leading to the same toroidal field 
strength through the fi-effect. Since the Lorentz force 
is proportional to the product of poloidal and toroidal 
field, the feedback on fl has almost doubled. This leaves 
however the question if it is possible to obtain larger 
toroidal field strength by going in the opposite direction 
and reducing the poloidal field strength. Inspecting Eq. 
© gives and estimate of the toroidal/poloidal field ratio, 
which can be obtained through the fl-eSect: 



t-^f--' 



(36) 



where we used r = 5 years and d^/dO ~ 0.2 J^o- For 
the model with ao = 0.125 ms~^ Bg within the convec- 
tion zone peaks at around 0.03 T, which does not allow 
for much larger toroidal field keeping in mind the reduc- 
tion of the poloidal field in this thought experiment and 
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considering that the estimate Eq. H36|l does not include 
resistive loss during the amplification phase. 

The cycle variation of fi (torsional oscillation) increases 
roughly proportional to ao from around 4.7 nHz to 18 
nHz at the pole, the change at 60° latitude is about half 
that value. Given the observational constraint that solar 
torsional oscillations at 60° have an amplitude of around 
2 nHz, the solar dynamo is most likely operating in a 
regime with weak Lorentz force feedback, close to our 
model with ag = 0.125 ms~^ or even smaller. In that 
model the amount of energy converted to magnetic en- 
ergy is around 8% of the energy converted into differen- 
tial rotation by the A-effect, in absolute units, around 
0.1% of the solar energy flux. This number is inter- 
estingly very close to the variation of solar irradiance 
throughout the solar cycle. 

In the last column of Table Q we show results ob- 
tained from a model with ag = 0.125 ms~^, but mag- 
netic diffusivity of only 5 x lO^m^s^^ in the convec- 
tion zone. This model is very close to the model with 
oq = 0.25ms~^ and a diffusivity of lO^m^s"^. The 
magnetic field strength as well as magnetic energy is 
slightly higher, at the same time the amount of energy 
extracted from differential rotation is lower due to the 
lower magnetic diffusivity. Nevertheless, the amplitude 
of torsional oscillations is larger compared to the model 
with higher diffusivity. 

3.3. Saturation field strength of dynamo 

Various models of rising magnetic fl ux - tubes 

dOiQudhuri fc GilnSp _ ITflSTt iFan et a,lJ [199^ 
ISchiissler et all Il994t ICaUgari et alJ 119951 119981) in- 
ferred a field strength of around lOT at the base of the 
convection zone, which is almost one magnitude more 
than the field strength we obtain in our model. Since our 
model is a mean field model, the mean field strength is 
not necessarily identical to the field strength of individ- 
ual flux elements in case of an intermittent fleld. When 
we compare our mean field solution to a solution with 
an intermittent field it is required to compare energy 
densities (which also ensure the same Lorentz force 
densities) rather than average field strength, since this 
ensures the same level of dynamic feedback. Suppose an 
intermittent field with a filling factor / and individual 
flux elements with field strength Bf. Conserving 
the mean energy leads a relation between mean field 
strength B and Bf of Bf ~ B/y/J. For a filling factor 
of / = 0.1 this would allow for around 4T for individual 
fiux elements. Another way to estimate maximum field 
strength is using the field energy directly. If we assume 
that the toroidal magnetic field is stored at the base of 
the convection zone close to the equator in a layer with 
a thickness Ar and a latitudinal extent Al, the magnetic 
energy is Eb ~ 2TTrArAlB^/{2fio) = nr^B/fiQ, with 
the magnetic flux $. Observations lead to an estimate 
of the ma gnetic flux of 10^ *^ - 10^^ Wb = 10^^ - 1 0^'^ M x 
llGallowav fc WeissI IT98lt ISchriiver fc HarvevI IT991 : 
however, the magnetic flux produced during a cycle 
at the base of the convection zone is not necessarily 
identical to the observed surface flux. On the one hand 
a rising flux rope could produce more than one spot 
group (which would allow for less flux at the base of the 
convection zone) on the other hand it is also likely that 
only a fraction of the flux at the base of the convection 



zone shows up in spots at the surface. The flux value 
we mention in Table ^ is the maximum flux of the 
dominant polarity available at a given time, while the 
estimates from observations are the cumulative values 
of flux integrated over the cycle. 

If we use the average energy conversion rate of 0.001 Fq 
and integrate that over a cycle length of 11 years we 
end up with an upper (very optimistic) estimate of 
1.3 X 10'^^ J, not considering any dissipative loss. In this 
case a value of <I> = 10^^ Wb would allow for i? = lOT, 
the higher value of $ = lO^'^Wb for only B = IT. 
Accounting for the additional dissipative loss, a field 
strength of more than a few T seems unlikely unless the 
average energy conversion rate would be larger. In our 
model that would lead to a contradiction with the ob- 
served amplitude of torsiona l oscillations. ^_^ 

It has been pointed out bv lRempel fc Schusslej l|200H) 
that also potential energy of the superadiabatic convec- 
tion zone can be used for the amplification of magnetic 
field. In contrast to the amplification through the fl- 
effect (shearing up of poloidal field by differential rota- 
tion), this does not lead to a feedback on differential 
rotation. Recently Y. Fan (2006, private communica- 
tion) repeated simulations of rising flux tubes using an 
anelastic 3D MHD code. These simulations are therefore 
not bound to the thin fiux tube approximation previ- 
ously used. Preliminary results show that the low lati- 
tude emergence and the observed asymmetries between 
leading and following spots (the tilt angle is a more com- 
plicated problem due to the infiuence of Coriolis force 
and twist of the fiux tube, which is currently being in- 
vestigated) can be reproduced with magnetic fiux tubes 
having an initial field strength only around 3 T, which 
would therefore relax the constraint on the field strength 
at the base of the convection zone. 



3.4. Role of magnetic buoyancy 

Our model includes magnetic buoyancy in two differ- 
ent ways: implicit buoyancy as part of the parametrized 
Babcock-Leighton a-effect and explicit (resolved) buoy- 
ancy resulting from solving the full momentum equation. 
While the first effect is essential for the dynamo model, 
the latter leads mainly to a distortion of the toroidal 
magnetic field shown in Figs. 01 to and is not essen- 
tial for the operation of the model. Having parametrized 
and resolved buoyancy together in a mean field model 
is of certain conceptual concern since it washes out the 
boundary between resolved and parametrized processes. 
Another point of concern is the fact that most buoyancy 
instabilities are non-axisymmetric and therefore our ax- 
isymmetric model does not capture the most unstable 
modes. One way of evaluating the importance of the 
explicit buoyancy in this model is to neglect the term 
'^ Pmag in Eq. Q that is responsible for it. This is 
equivalent to adding to si a perturbation of — Pmag/po, 
which means physically that the toroidal field is stored in 
a buoyancy free equilibrium at the base of the convection 
zone. The Babcock-Leighton a-effect addresses buoy- 
ancy instabilities below the resolved scale and is there- 
fore not in contradiction with this assumption. Since 
this additional entropy perturbation is not considered in 
Eq. (|SJ) this leads to a small inconsistency in the entropy 
equation. For the analysis of the energy fluxes we use 
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Fig. 8. — Magnetic field evolution for models with magnetic 
buoyancy switclicd off. Similar to Fig. |n] solutions with ag = 
0.125ms~^ (top), 0.25ms~^ (middle), and0.5ms~^ (bottom) are 
shown. 



instead of Eq. (|^ 



Qb = - dVvrQog- Si ^^^ 1 

J 7 V Po / 



(37) 



Fig. |H1 shows the magnetic field evolution similar to Fig. 
(0 with magnetic buoyancy switched off. While the top 
panel (solution with ao = 0.125 ms^^) is very close to 
the kinematic reference solution with a-quenching, the 
middle and bottom panels show dynamo solutions with 
irregular cycles due to the non-linear feedback. Despite 
the significant difference in the magnetic field pattern 
and temporal evolution, the differences in the energy ex- 
change terms are not that significant. The most obvious 
differences occur in the model with ao = 0.5 ms~^ in 
terms of amplitude of torsional oscillations, radial surface 
field strength, and variability of magnetic energy. The 
almost a factor of 2 larger radial field is also the reason 
for the large irregularity of the solution. The stronger 
poloidal field at the surface leads through the Lorentz 
force to a stronger variation in the meridional flow that 
changes the latitude at which the magnetic field is trans- 
ported downward. This in return changes the latitudinal 
extent and strength of the next cycle in a way that a 
periodic solution is not possible anymore. If buoyancy 




Fig. 9. — Influence of meridional flow profile on phase relation 
between high latitude torsional oscillation and low latitude toroidal 
magnetic field. Top: solution shown in Fig. |^ (bottom panel) for 
reference; bottom: model with meridional flow returning in higher 
latitudes (n=2 instead of n=3 in reference model). 



is considered the radial surface field does not reach the 
threshold required for a highly irregular solution. 

4. TORSIONAL OSCILLATIONS 

Solar torsional oscillations have been known to exist 
for more than two decades. iHoward fc Labon tc (198(4) 
presented the first observations of torsional oscillations 
using Mt. Wilson Doppler measurements and pointed 
out the 11 year periodicity and the relation to the so- 
lar cycle. These early observations showed only the 
equatorward propagating branch at low latitudes. The 
high latitude branch (above 60°), which is in amplitude 
at least twice as strong as the equatorward propagat- 
ing branch, was foun d more recently through helioseis- 
mic me asurements b v^Toomre e t al.l (12 0001): iHowe et a,l] 
(200d): lAntia fc Ba su (2001); Vorontsov et alJ (|2002|) : 
iHowe et al.l ()200 5l)7 These inversions also show that the 
high latitude signal penetrates almost all the way to the 
base of the convection zone. The depth penetration of 
the low latitude signal is more uncertain due to the lower 
amplitude that is comparable to the uncertainties of the 
inversion methods in the lower half of the convection 



4.1. Mechanical forcing of torsional oscillations 

The models discussed so far can only explain the po- 
lar branch of the torsional oscillation pattern as a very 
robust result through the feedback of mean field Lorentz 
force on differential rotation (mechanical forcing). We 
also showed that some dynamo models lead to quite sig- 
nificant amplitudes of these oscillations, contrary to ob- 
servation, which therefore impose constraints on dynamo 
parameters (amplitude of a-effect and the magnetic dif- 
fusivity in convection zone) . Other information that can 
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be extracted from torsional oscillations is their phase rel- 
ative to t he magnetic cycl e . The phase relation shown in 
Fig. 3 of iVorontsov et all l|2002f) is such that the pole is 
rotating faster during solar minimum and slower during 
maximum, when solar activity peaks around 15° latitude. 
In order to determine the phase relation in our model, 
we have to define what phase corresponds to 'solar min- 
imum' and 'solar maximum', given a butterfly diagram 
which is close to, but not exactly, solar like. If we define 
'solar maximum' through the maximum field strength at 
the base of the convection zone, the location of the ac- 
tivity belt is at around 40° and the torsional oscillation 
pattern changes from faster to slower rotation at the pole 
during that time. If we however define 'solar maximum' 
as the time where the activity belt is around 15° latitude, 
the phase relation of the model matches the observed 
torsional oscillations. Since the phase in a flux-transport 
dynamo is a consequence of the advection of field by the 
meridional flow, while the amplitude depends on details 
of microphysics that might not be well represented or 
even missing in this mean field model (e.g. a tachocline 
a-effect), the latter definition is most likely more robust 
and should be used when comparing results to the sun. 
The phase relation is indeed tied to the structure of the 
meridional flow. Fig. compares the model shown in 
Fig. 01 with model that has a meridional flow returning 
at higher latitudes (we used instead of n = 3 and Aq = 1 
n = 2 and Aq = 0.8 in the reference model). While in 
the original model slower rotation at the pole coincided 
with an activity belt location at around 20°, it coincides 
with an activity belt location at around 30° in the latter 
model. Also note that the dynamo period of the latter 
model is around 22.5 years instead of 18 years, due to 
the longer overturning time of the meridional flow. 

4.2. Thermal forcing of torsional oscillations 

An alternative explanation for the low l atitude branch 
of the torsional oscillations was given by Sorviit l'i 2003t) : 
Enhanced surface cooling in the active region belt leads 
to a pressure imbalance that drives a geostrophic flow at 
the edges of the active region belt showing the properties 
of the observed pattern . The low latitude oscillation pat- 
tern would be therefore purely surface driven and not a 
signature of Lorentz force feedback within the convection 
zone. 

We can include the idea of Spruit in our model by 
parameterizing a surface cooling term that depends on 
the toroidal field at the base of the convection zone. To 
this end we change the upper boundary condition such 
that the entropy gradient corresponds to an increase in 
energy flux. With the diffusive convective energy flux 
Fc = ntgoToCyVsi and c„ = (7 - l)"Vo/(£'o To) this 
leads to 



9s 1 
dr 



7- 



1 



KtPO 



■eFr, 



(sin6')^.B^^bc 

-Brcf 



(38) 



Here i?$,bc(^) is the toroidal field averaged according to 
Eq. (fl^ between 0.71 i?© and 0.76 i?©, which is also used 
for the Babcock-Leighton a-effect. B^^i is a reference 
field strength used for normalization, e determines the 
amplitude of the flux enhancement as fraction of the so- 
lar luminosity. In order to illustrate this effect and allow 
a comparison to the observed solar torsional oscillation 




Fig. 10. — Torsional oscillations caused through enhanced ra- 
diative losses in the active region belt. Top: reference model with 
no surface cooling; middle: model with surface cooling, the maxi- 
mum amplitude is 4nHz, bottom: surface temperature variation, 
the maximum amplitude is 0.2 K. 



we introduce a factor (sin 6)"^ in front of B^ to correct the 
latitudinal field strength profile of the butterfiy diagram 
accordingly (this is consistent with the {sin 9)^ factor in- 
troduced in the a-effect for the same reason). We use a 
dynamo model with magnetic buoyancy switched off and 
a smaller value of a (0.1 ms~^) than before in order to 
obtain an amplitude of the torsional oscillations close to 
the observed one (around 4 nHz at the pole and 2 nHz at 
60° latitude. Fig. ^| (top panel) shows the torsional os- 
cillation and butterfly diagram for the dynamo solution 
not considering the surface cooling. The middle panel 
shows the solution with surface cooling considered using 
the parameters e = 2 x 10^^ and B^^f = 1 T. The cooling 
of the active region belt drives a equatorward propagat- 
ing torsional oscillation with an amplitude of around 1.5 
nHz, with the peak values at the edges of the active re- 
gion belt. The amplitude of the associated temperature 
perturbations (shown in the bottom panel) is around 0.2 
K. The peak cooling rate in our model is close to 0.75% of 
the solar energy flux, the surface intergrated luminosity 
variation corresponds to around 0.23% luminosity change 
throughout the cycle, about a factor of 3 larger than the 
observed one. We emphasize that this value is imposed 
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Fig. 11. — Torsional oscillations for model with additional surface 
cooling. The frames (top to bott om) correspond to the times t = 
1, t = 3.5, and t = 6, in Fig. 1111 While the polar branch of 
the torsional oscillations penetrates all the way to the base of the 
convection zone, the equatorial branch caused by surface cooling is 
confined close to the surface. 
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Fig. 12. — Top: near surface torsional oscillation caused through 
quenching of turbulent viscosity, the maximum amplitude is 0.77 
nHz; bottom: torsional oscillation caused through quenching of 
turbulent viscosity and heat conductivity, the maximum amplitude 
is 3.1 nHz, respectively. 



at r = 0.985 Rq and does not necessarily resemble the 
value required in a more realistic model extending all the 
way into the photosphere and using a more sophisticated 
description of convection than the diffusion approxima- 
tion. Fig. E] shows three snapshots of the evolution of 
the torsional oscillation in a r-0-plane. While the pole- 
ward propagating branch (driven by Lorentz force feed- 
back on differential rotation) is penetrating all the way 
down to the base of the convection zone, the equatorward 
propagating branch (driven through the surface cooling) 
is more concentrated toward the surface. As a side effect 
of the surface cooling the model also shows close to the 
surface an inflow into the active region belt with a peak 
amplitude of aroun d 5 ms ~ ^. Th is value is in agreement 
with the theory of lSpruilJ l|2003j) . where the meridional 
component is a conse quence of an Ekman boundary layer 
at the solar surface. iKomm et"aLl l|1993() : iKomml l)1994|) 
derived an average inflow into the active region belt of 
around 5ms^^ from Kitt Peak magnetograms. More re- 
cently also helioseisniology showed a mean inflow into 
the active region belt with an amplitude from 2ms^^ to 



8ms-i (|Zhao fc Kosovichevll20& 



4.3. Forcing of torsional oscillations through quenching 
of viscosity and heat conductivity 

So far we only considered the 'macroscopic' Lorentz 
force feedback in terms of the Lorentz force computed 
from the magnetic mean fleld. Alternatively we can con- 
sider 'microscopic' feedback through quenching of tur- 
bulent motions. This type of feedback is typically used 
in kinematic models in terms of a-quenching to saturate 
the dynamo. We consider here the quenching of turbu- 
lent viscosity and heat conductivity. Since in our model 
the turbulent viscosity scales the amplitude of the A- 
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effect, quencliing of vt reduces the energy input into the 
system and leads to a saturation of the dynamo as con- 
sequence of a reduction in differential rotation similar 
to the 'macroscopic' feedback discussed above. Addi- 
tional quenching of turbulent heat conductivity leads to 
changes in entropy profile of the convection zone. Since 
the differential rotation is close to a baroclinic balance, 
a change in the entropy profile also forces changes in the 
differential rotation. Fig. ^| (top panel) shows torsional 
oscillations and the butterfiy diagram obtained from a 
dynamo model with ag = 0.125 ms~^ and Vt quenching 
with Beq = 1 T according to Eq. JTHl. Fig. [l3 (bottom 
panel) shows torsional oscillations resulting from addi- 
tional quenching of the turbulent heat conductivity. Note 
the different phase relation of the torsional oscillations 
with respect to the magnetic butterfly diagram compared 
to the results with 'macroscopic' Lorentz force feedback. 
Quenching of turbulent heat conductivity alone is not a 
process that can efficiently saturate the dynamo, since 
it leads more to a modulation of fJ rather than a reduc- 
tion in the equator-pole difference. It can however sig- 
nificantly change the torsional oscillation pattern. Since 
there are additional effects such as anisotropic heat trans- 
port, which are not considered in our model, we empha- 
size here only the interesting result that the amplitude 
of these thermally forced oscillations is quite significant. 
Quenching of turbulent viscosity is along the lines of 
A-que n ching considered b efore b v iKitchatinov fc Pioinl 
l|T998|) : lKitcliatinov et all IpOQliK fiker et all l|19m 

5. PARAMETER DEPENDENCE 

The results presented in this paper were obtained us- 
ing a mean field model that requires parameterizations of 
unresolved processes. Therefore it is necessary to test to 
which extent results are sensitive with respect to details 
of the parameterizations used. Our differential rotation / 
meridional fiow model has three important parameters: 
n describing the profile of the A-effect in latitude, Aq 
determining the amplitude of the A-effect, and the tur- 
bulent diffusivities vt, f^t (both have same value in our 
model). Aq has been chosen to get differential rotation 
with the correct amplitude. For a fixed value of Aq, vt 
and Kt determine the amplitude of the meridional flow. 
The latter has been chosen to lead to dynamo simula- 
tions with a period close to 22 years. The parameter n 
determines the extent of the meridional flow cell in lati- 
tude. Reasonable choices are in the range 2—4 (for values 
above 4 differential rotation is confined to low latitudes, 
which also contradicts observations). Within this range 
the influence on the solutions is rather limited (see Fig. 
in|) . A fourth parameter not discussed in this paper is the 
direction of the turbulent angular momentum flux with 
respect to the axis of rotation. In this paper we used a 
fixed valu e of A = 15°, as has been used also for most 
models in iRemoell l)2nn 5b') . Significantly larger values 
lead to more complicated meridional fiow patterns (per- 
manent reverse cell in high latitudes), smaller values re- 
quire and increase of Ag to maintain the differential rota- 
tion amplitude, which in return leads to larger meridional 
flow speeds, requiring signiflcantly lower values for i^t and 
Kt to obtain a solar-like dynamo period. We computed 
solutions (not shown in this paper) with the parame- 
ters A = 7.5°, An = 2 and vt = Kt = 1.25 x lO^m^s^^ 
that show no significant difference, except that the am- 



plitude of the torsional oscillations is around 50% larger 
due to the reduced viscous damping. Therefore calibrat- 
ing our reference model to be solar-like does not leave a 
lot of choice for the parameters of the differential rota- 
tion model. Using such a calibrated reference model, the 
following main results are not strongly dependent on the 
additional dynamo parameters (a-effect, rjt) providing 
the magnetic diffusivity is low enough to be in an ad- 
vection dominated regime: The dynamo saturates with 
a toroidal (mean) field strength of about 1.2 to 1.5 T at 
the base of the convection zone, the maximum magnetic 
flux produced during a cycle is around 10^^ Wb, the en- 
ergy converted on average into magnetic energy is 0.1% 
to 0.2% of the solar luminosity (more likely 0.1 incorpo- 
rating additional constraints set by helioseismic observa- 
tions of torsional oscillations). Another important result 
is that the dynamo saturates already through reduction 
of differential rotation at a field strength where the direct 
feedback on the meridional flow is insigniflcant. There- 
fore the fundamental character of the flux-transport dy- 
namo remains unaltered (the non-linear feedback tends 
to concentrate magnetic activity even closer toward the 
equator). 

6. IMPLICATIONS FOR SOLAR DYNAMO 

In this paper we discussed 'dynamic' flux-transport dy- 
namos by combining a mean field differential rotation 
model with a mean field dynamo model. The main re- 
sults from this study are the following 

• The non-kinematic Babcock-Lcighton flux- 
transport dynamo saturates through a reduction 
of the amount of differential rotation at a toroidal 
field strength of around 1.5 T. This field strength 
is found to be fairly independent from particular 
choices of the strength of the a-effect and value of 
magnetic diffusivity. The energy conversion rate 
of the dynamo is around 0.001 Lq. 

• The Lorentz force feedback (at the saturation field 
strength of around 1.5 T) on the meridional flow 
is not strong enough to switch off the equatorward 
transport of toroidal fleld required for the operation 
of the flux-transport dynamo . This is consistent 
with the findings of iRemnell l|200fil) who studied 
the modiflcation of meridional flow by an imposed 
stationary toroidal magnetic field. 

• The Lorentz force feedback on differential rotation 
leads to torsional oscillations that are comparable 
(in terms of amplitude and phase with respect to 
the magnetic butterfiy) to the high latitude branch 
of the pattern inferred from helioseismology. The 
low latitude branch cannot be explained through 
Lorentz force feedback in our model. A model in- 
cluding thermal forcing through increased radiative 
losses in the active region belt (.Spruit :2,Q0.3) also 
produces the low latitude branch. Thermal forcing 
is found to be very efficient in the sense that tem- 
perature variation of only a few tenth of a degree 
are sufficient to explain the observed amplitudes of 
torsional oscillations. 

The results presented here cannot answer the ques- 
tion of whether the solar dynamo is a flux-transport dy- 
namo or not; however, they make the case stronger for a 
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flux transport dynamo since they demonstrate that flux- 
transport dynamos also function in the non-kinematic 
regime (the meridional flow is strong enough to trans- 
port toroidal field of around 1.5 T strength equatorward) . 
Additional to that the torsional oscillations caused by 
the macroscopic Lorentz force in a flux-transport dy- 
namo are in agreement with helioseismic results (in terms 
of amplitude and phase relation). With 'in agreement' 
we mean here that there are no features produced that 
are not observed; however there are observed features 
(the low latitude branch ) that r equire additional physics 
such as thermal forcing IjSDruitl [20031. The idea of ther- 
mal forcing of the low latitude branch is also consistent 
with a mean inflow into the active region belt of the or- 
der of 5ms~^ that has been observed in magnetograms 
IjKomm et al.lU,993: Ko mm 1994) and in ferred from he- 
lioseismology llZhao fc K osovichevl l200^ . We emphasize 
that thermal forcing is a very efficient process, since only 
tiny temperature fluctuations of the order of a few tenth 
of a degree can drive large scale flows with the observed 
amplitude. This opens the possibility that also other pro- 
cesses such as the magnetic quenching of the convective 
e nergy flux in the solar convection zone contribute. 

ICovaset al1 l2000. 2004, 2005) presented a al7-dynamo 
model (no meridional flow) including a simplified mo- 
mentum equation for considering feedback on differential 
rotation. In their simulation they were able to repro- 
duce additional to the polar branch also the low latitude 
oscillations pattern. The main difference between their 
and our model is that we use a flux-transport dynamo 
with a non-local Babcock-Leighton a-effect, while they 
use a classic arj-dynamo model with a (negative) local 
Qf-effect in the convection zone to obtain an equatorward 
propagating dynamo wave. The requirement to have a 
propagating dynamo wave leads to a fixed phase relation 
between poloidal and toroidal field, which automatically 
leads to a Lorentz force pattern propagating with the 
field and producing a torsional oscillation pattern asso- 
ciated with the magnetic field. This constraint is relaxed 
in a flux-transport dynamo, where the propagation of 
the activity belt is a pure advection effect. Given the 
fact that the origin of the low latitude oscillation pat- 
tern in uncertain and might be entirely surface driven, it 
is difficult to judge to which extent this discrepancy is of 
concern for different types of dynamo models and might 
rule out certain approaches. 

We emphasize here that it is essential for models of 
torsional oscillations to include the full momentum equa- 
tion, since the Taylor-Proudman constraint applies also 
to perturbations of fi {di^i/dz = 0) and therefore signif- 
icantly alters the phase relation of the oscillations. Tor- 
sional oscillations with dVti/dz ^ (as observed in the 
sun at low latitudes) require additional thermal pertur- 
bations, which makes also the consideration of an energy 



equation essenti al and f avors explanations such as the 
one proposed bv ISorui^ lj2003|) . 

Torsional oscillations contain valuable information 
about the dynamo processes in the solar interior and 
are helpful to impose additional constraints on dynamo 
models. Their interpretation and relation to the dynamo 
generated magnetic field is complicated since different 
processes can produce torsional oscillations. As a first 
step toward using torsional oscillations to probe solar 
cycle related processes in the solar interior it is impor- 
tant to test the sensitivity of helioseismic inversions to 
dynamo generated rotation modulations in the solar in- 
terior, especially close to the ba se of the con vection zone. 
This work is currently done bv lHowe et alJ (|^04, 200f| 
for different helioseismic techniques using artificial data 
and also model results discussed in this paper. 

The model presented here sets strong constraints for 
the strength of toroidal field at the base of the convec- 
tion zone that can be achieved through the shear by 
differential rotation. We find in our model a value of 
around 1.5 T (15 kG) as upper limit. These resu lts are 
in agr eement with recent findings of Gilman & R empel 
(J2005J) who showed that significantly larger field strength 
would require a very strong mechanism replenishing en- 
ergy to differential rotation (replenishment time-scale of 
less than a month). The fairly low field strength (com- 
pared to convective equipartition) raises the question of 
how these fields can rise through the convection zone 
and form coherent sunspots at the surface. On the one 
hand recent work by Y. Fan (2006, private communi- 
cation) suggests that 3D simulations of rising flux tubes 
can reproduce most sunspot properties starting with field 
strengths around 2 to 3 T at the base of the convection 
zone (as opposed to thin flux tube simulations requir- 
ing arovm d 10 T ); on the other hand recent work by 
Brummcll et al.l (|2002); Cattanco et al. (2006) is ques- 
tioning the existence of coherent buoyant structures. In 
our model these uncertainties are hidden behind the pa- 
rameterization of the Babcock-Leighton a-effect. Since 
the existence of this effect is known from surface observa- 
tions (it is the essential ingredient in models describing 
the evolution of th e surface magnetic field assimilating 
real observ ations (.Schriiver et 311120021 IBaumann et alJ 
I2OOI I20M IWang et al.ll2005fl ) this is of secondary con- 
cern for the model presented here; it is however of pri- 
mary interest for understanding the microphysics beyond 
the mean field approach. 
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